# -*- coding: utf-8 -*-
"""
:Author: Fernando Pacheco <fernando.pacheco@ingesur.com.uy>
:Date: Fri, 10 Jul 2009 12:31:51 -0300
:copyright: Ingesur SRL
:license: GPL Version 2
"""
import math
from scipy import integrate
from scipy import special

class BaseFunctions(object):
    """
    Base function class. Common function in hydrogeology.
    
    :note: Some functions in the class have very particular names (without under case)
    """

    def __init__(self):
        pass

    def _fu(self, x):
        """
        Integration function for Wu (confinated and unconfinated aquifers)
            Wu=Integral(fu) between u and +Inf
        
        :Parameters:
            x : float 
                Integral variable

        """
        return (1/x*math.exp(-x))

    def wu_integrate(self, u):
        """
        Integral of function fu between u and ls.
        Calculated with integrate function of numpy package.
        Usually used in confined and unconfined aquifers in nonpermantent regime.

        :Parameters:
            u : float
                u is r^2*S/(4*T*t)
        """
        return integrate.quad(self._fu,u,integrate.Inf)

    def wu_series(self, u):
        """
        Series expansion of Wu function.
        Usually used in confined and unconfined aquifers in nonpermantent regime.
        
        :Parameters:
            u : float
                u is r^2*S/(4*T*t)
        """
        if u <= 0: 
            return ( [0, 0] )            
        eps = 0.0001
        su = - 0.577215665 - self.log(u)
        if u <= eps:
            return ( [su, u] )            
        tt = (-1)*u
        su = su - tt            
        if u < 1:
            n = 2
            while abs( tt ) >= eps:
                tt = (-1)*tt*u*(n-1)/(n*n)
                su = su - tt 
                n = n+1
                # print "n=%d,u=%f,tt=%f,su=%f" % (n, u, tt, su) 
            return ( [su, u] )    
        else:
            if u<=21:
                # Con N ~ 4*(u+1) la convergencia es bastante buena
                for cont in range(2, 4*(u+1)):
                    tt = (-1)*tt*u*(cont-1)/( float(cont*cont) )
                    su = su - tt            
                    # print "cont=%d,u=%E,tt=%E,su=%E" % (cont, u, tt, su) 
                return ( [su, u] )
            else:
                # Arroja valores muy bajos y es demasiado costoso calcularlo
                # Ej. u = 18 Wu = 6.768184E-10
                return ( [eps*0.00001, u] )

    def _fua(self, x, a):
        """
        Integration function for Wua (semiconfined aquifers)  
        Wua=Integral(fua) between u and +Inf
    
        :Parameters:
            x : float
                Integral variable
        """
        return (1/x*math.exp(-x-a*a/(4*x)))
    
    def wua_integrate(self, u, a):
        """
        Function for calculate fua integral. Well function for semiconfined aquifers.
        Calculated with integrate function of numpy package.
    
        :Parameters:
            u : float
                Donde u es r^2*S/(4*T*t)
            a : float
                (alpha) Donde alpha es r/B
        """
        return integrate.quad(self._fua,u,integrate.Inf, args = (a))        

    def k0(self, x):
        """
        Bessel function 0 order. Calculated with k0 function of numpy package.
        Usually used in semiconfined aquifers (permanent regime) 
        
        :Parameters:
            x : float
                Variable (x = r/B)
                
        """        
        return(special.k0(x))
        
    def k1(self, x):
        """
        Bessel function 1 order. Calculated with k1 function of numpy package.
        Usually used in semiconfined aquifers (permanent regime) 
        
        :Parameters:
            x : float
                Variable (x = r/B)                
        """        
        return(special.k1(x))
        
    def pi(self):
        """
        Pi number
        """
        return(math.pi)

    def log(self, x):
        """
        Natural logarithmic, log(10) = 2.3025850929940459
        """
        return(math.log(x))
        
    def exp(self, x):
        """
        Natural logarithmic, e(1) = 2.718281828459045
        """
        return(math.exp(x))
        
    def sqrt(self, x):
        """
        sqrt(4) = 2
        """    
        return(math.sqrt(x))

    def linear_interpolation(self,x,pa,pb):
        """
        Linear interpolation.

        :Parameters:
            x : float
                Abscisa value to calculate
            pa : list of float
                (xa,ya) = values
            pb : list of float
                (xb,yb) = values
        """
        if (x > pa[0] and x > pb[0]) or (x < pa[0] and x < pb[0]):
            print "Interpolation error."
            val = None
        else:
            val = pa[1] + ((x - pa[0]) * (pb[1] - pa[1]))/(pb[0] - pa[0])
        return val